comparative_analysis_TGFbhigh_vs_low

workflowr

• Summary
• Checks
• Past versions

Last updated: 2024-08-27

Checks: 5 2

Knit directory: CART_TGFb/

This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

---

R Markdown file: uncommitted changes

The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20240213)

The command set.seed(20240213) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: absolute

Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.

absolute	relative
/home/hnatri/CART_TGFb/code/CART_plot_functions.R	code/CART_plot_functions.R
/home/hnatri/CART_TGFb/code/13384_tumor_ms_themes.R	code/13384_tumor_ms_themes.R
/home/hnatri/CART_TGFb/code/utilities.R	code/utilities.R
/home/hnatri/CART_TGFb/	.
/home/hnatri/CART_TGFb/data/ELISA1_TGFb.tsv	data/ELISA1_TGFb.tsv

Repository version: 33685f1

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 33685f1. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:

Ignored files:
    Ignored:    analysis/figure/
    Ignored:    code/.RData

Unstaged changes:
    Modified:   analysis/analyze_compartments.Rmd
    Modified:   analysis/comparative_analysis_TGFbhigh_vs_low.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

---

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/comparative_analysis_TGFbhigh_vs_low.Rmd) and HTML (docs/comparative_analysis_TGFbhigh_vs_low.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	33685f1	heinin	2024-08-27	Moving files
Rmd	9015511	heinin	2024-08-26	comparative analysis and TGFb expression in scRNAseq
html	9015511	heinin	2024-08-26	comparative analysis and TGFb expression in scRNAseq

---

Comparative analysis of immune lineages between TGFb high/low

Libraries, helper functions, and environment variables

library(workflowr)
library(Seurat)
library(googlesheets4)
library(dplyr)
library(tidyverse)
library(ggrepel)
library(patchwork)
library(scProportionTest)
library(escape)

source("/home/hnatri/CART_TGFb/code/CART_plot_functions.R")
source("/home/hnatri/CART_TGFb/code/13384_tumor_ms_themes.R")
source("/home/hnatri/CART_TGFb/code/utilities.R")

setwd("/home/hnatri/CART_TGFb/")

set.seed(1234)
options(future.globals.maxSize = 30000 * 1024^2)
reduction <- "umap"

Import data

immune_seurat <- readRDS("/tgen_labs/banovich/BCTCSF/Heini/13384_Tumor/immune_reclustered_TGFb.rds")
tumor_seurat <- readRDS("/tgen_labs/banovich/BCTCSF/Heini/13384_Tumor/tumor_reclustered_TGFb.rds")

# Adding skull cluster numbers
immune_robyn <- readRDS("/scratch/hnatri/CART/fromRobyn/IMMUNEminus618_TGFbELISA+TGFbELISA2.rds")

DimPlot(immune_robyn, group.by = "seurat_clusters", label = T)

Past versions of unnamed-chunk-2-1.png

Version	Author	Date
9015511	heinin	2024-08-26

length(intersect(colnames(immune_seurat), colnames(immune_robyn)))

[1] 25098

length(setdiff(colnames(immune_seurat), colnames(immune_robyn)))

[1] 504

head(immune_robyn@meta.data)

                             orig.ident nCount_RNA nFeature_RNA percent.mt
Batch1_ACGCAGCAGTCATCCA-1 SeuratProject      10686         2552   2.227213
Batch1_ACTGAGTGTTCCGGCA-1 SeuratProject       4613         1652   2.124431
Batch1_AGCATACGTAGCTAAA-1 SeuratProject       7305         2086   3.449692
Batch1_ATGTGTGCAATAGCGG-1 SeuratProject       3776         1533   4.555085
Batch1_CAGTCCTAGCTTTGGT-1 SeuratProject       2719         1457   9.157779
Batch1_CGATTGATCTATCCCG-1 SeuratProject       9146         2694   4.089219
                          nCount_Protein nFeature_Protein
Batch1_ACGCAGCAGTCATCCA-1           4689              188
Batch1_ACTGAGTGTTCCGGCA-1           3599              187
Batch1_AGCATACGTAGCTAAA-1           3095              179
Batch1_ATGTGTGCAATAGCGG-1           1557              180
Batch1_CAGTCCTAGCTTTGGT-1           4541              191
Batch1_CGATTGATCTATCCCG-1           4598              191
                                Exome_Sample_Name Demultiplex_Assignment UPN
Batch1_ACGCAGCAGTCATCCA-1 BCTCSF_0058_1_PB_MNC_C1                    SNG 243
Batch1_ACTGAGTGTTCCGGCA-1 BCTCSF_0058_1_PB_MNC_C1                    SNG 243
Batch1_AGCATACGTAGCTAAA-1 BCTCSF_0058_1_PB_MNC_C1                    SNG 243
Batch1_ATGTGTGCAATAGCGG-1 BCTCSF_0058_1_PB_MNC_C1                    SNG 243
Batch1_CAGTCCTAGCTTTGGT-1 BCTCSF_0058_1_PB_MNC_C1                    SNG 243
Batch1_CGATTGATCTATCCCG-1 BCTCSF_0058_1_PB_MNC_C1                    SNG 243
                          Sample_Type Cycle Day Manufacture     Cycle_Day
Batch1_ACGCAGCAGTCATCCA-1         FDT    NA  NA       Tumor CycleNA_DayNA
Batch1_ACTGAGTGTTCCGGCA-1         FDT    NA  NA       Tumor CycleNA_DayNA
Batch1_AGCATACGTAGCTAAA-1         FDT    NA  NA       Tumor CycleNA_DayNA
Batch1_ATGTGTGCAATAGCGG-1         FDT    NA  NA       Tumor CycleNA_DayNA
Batch1_CAGTCCTAGCTTTGGT-1         FDT    NA  NA       Tumor CycleNA_DayNA
Batch1_CGATTGATCTATCCCG-1         FDT    NA  NA       Tumor CycleNA_DayNA
                           Batch nCount_SCT nFeature_SCT     S.Score
Batch1_ACGCAGCAGTCATCCA-1 Batch1       3133         1091  0.02190489
Batch1_ACTGAGTGTTCCGGCA-1 Batch1       3516         1647 -0.01952535
Batch1_AGCATACGTAGCTAAA-1 Batch1       3639         1679 -0.04382463
Batch1_ATGTGTGCAATAGCGG-1 Batch1       3385         1533 -0.02855645
Batch1_CAGTCCTAGCTTTGGT-1 Batch1       2821         1457 -0.09747976
Batch1_CGATTGATCTATCCCG-1 Batch1       3526         1709 -0.06602777
                             G2M.Score Phase   IRB  Sample_ID nCount_SoupX_RNA
Batch1_ACGCAGCAGTCATCCA-1  0.008389157     S 13384 243_Batch1                0
Batch1_ACTGAGTGTTCCGGCA-1 -0.031608446    G1 13384 243_Batch1                0
Batch1_AGCATACGTAGCTAAA-1  0.036195260   G2M 13384 243_Batch1                0
Batch1_ATGTGTGCAATAGCGG-1 -0.108682135    G1 13384 243_Batch1                0
Batch1_CAGTCCTAGCTTTGGT-1  0.045369327   G2M 13384 243_Batch1                0
Batch1_CGATTGATCTATCCCG-1 -0.012711941    G1 13384 243_Batch1                0
                          nFeature_SoupX_RNA percent.mt_RNA percent.ribo_RNA
Batch1_ACGCAGCAGTCATCCA-1                  0             NA               NA
Batch1_ACTGAGTGTTCCGGCA-1                  0             NA               NA
Batch1_AGCATACGTAGCTAAA-1                  0             NA               NA
Batch1_ATGTGTGCAATAGCGG-1                  0             NA               NA
Batch1_CAGTCCTAGCTTTGGT-1                  0             NA               NA
Batch1_CGATTGATCTATCCCG-1                  0             NA               NA
                          percent.mt_SoupX_RNA percent.ribo_SoupX_RNA
Batch1_ACGCAGCAGTCATCCA-1                   NA                     NA
Batch1_ACTGAGTGTTCCGGCA-1                   NA                     NA
Batch1_AGCATACGTAGCTAAA-1                   NA                     NA
Batch1_ATGTGTGCAATAGCGG-1                   NA                     NA
Batch1_CAGTCCTAGCTTTGGT-1                   NA                     NA
Batch1_CGATTGATCTATCCCG-1                   NA                     NA
                          nCount_Hash nFeature_Hash Hash_maxID Hash_secondID
Batch1_ACGCAGCAGTCATCCA-1           0             0       <NA>          <NA>
Batch1_ACTGAGTGTTCCGGCA-1           0             0       <NA>          <NA>
Batch1_AGCATACGTAGCTAAA-1           0             0       <NA>          <NA>
Batch1_ATGTGTGCAATAGCGG-1           0             0       <NA>          <NA>
Batch1_CAGTCCTAGCTTTGGT-1           0             0       <NA>          <NA>
Batch1_CGATTGATCTATCCCG-1           0             0       <NA>          <NA>
                          Hash_margin Hash_classification
Batch1_ACGCAGCAGTCATCCA-1          NA                <NA>
Batch1_ACTGAGTGTTCCGGCA-1          NA                <NA>
Batch1_AGCATACGTAGCTAAA-1          NA                <NA>
Batch1_ATGTGTGCAATAGCGG-1          NA                <NA>
Batch1_CAGTCCTAGCTTTGGT-1          NA                <NA>
Batch1_CGATTGATCTATCCCG-1          NA                <NA>
                          Hash_classification.global hash.ID FID_GEXFB
Batch1_ACGCAGCAGTCATCCA-1                       <NA>    <NA>      <NA>
Batch1_ACTGAGTGTTCCGGCA-1                       <NA>    <NA>      <NA>
Batch1_AGCATACGTAGCTAAA-1                       <NA>    <NA>      <NA>
Batch1_ATGTGTGCAATAGCGG-1                       <NA>    <NA>      <NA>
Batch1_CAGTCCTAGCTTTGGT-1                       <NA>    <NA>      <NA>
Batch1_CGATTGATCTATCCCG-1                       <NA>    <NA>      <NA>
                          seurat_clusters   IFNg_score_SER IFNg_score_TF
Batch1_ACGCAGCAGTCATCCA-1              10 2.64759961886982          <NA>
Batch1_ACTGAGTGTTCCGGCA-1               1 2.64759961886982          <NA>
Batch1_AGCATACGTAGCTAAA-1               5 2.64759961886982          <NA>
Batch1_ATGTGTGCAATAGCGG-1               9 2.64759961886982          <NA>
Batch1_CAGTCCTAGCTTTGGT-1               4 2.64759961886982          <NA>
Batch1_CGATTGATCTATCCCG-1              10 2.64759961886982          <NA>
                          IFNg_score_CSF Survival.time.in.Months.from.surgery
Batch1_ACGCAGCAGTCATCCA-1           <NA>                                   NA
Batch1_ACTGAGTGTTCCGGCA-1           <NA>                                   NA
Batch1_AGCATACGTAGCTAAA-1           <NA>                                   NA
Batch1_ATGTGTGCAATAGCGG-1           <NA>                                   NA
Batch1_CAGTCCTAGCTTTGGT-1           <NA>                                   NA
Batch1_CGATTGATCTATCCCG-1           <NA>                                   NA
                          Death.Status DS Grade Diagnosis.Histology
Batch1_ACGCAGCAGTCATCCA-1           NA NA  <NA>   Glioblastoma, NOS
Batch1_ACTGAGTGTTCCGGCA-1           NA NA  <NA>   Glioblastoma, NOS
Batch1_AGCATACGTAGCTAAA-1           NA NA  <NA>   Glioblastoma, NOS
Batch1_ATGTGTGCAATAGCGG-1           NA NA  <NA>   Glioblastoma, NOS
Batch1_CAGTCCTAGCTTTGGT-1           NA NA  <NA>   Glioblastoma, NOS
Batch1_CGATTGATCTATCCCG-1           NA NA  <NA>   Glioblastoma, NOS
                          PriorAvastin No.Prior.Chemo No.Prior.Brain.Surgeries
Batch1_ACGCAGCAGTCATCCA-1          yes              3                        2
Batch1_ACTGAGTGTTCCGGCA-1          yes              3                        2
Batch1_AGCATACGTAGCTAAA-1          yes              3                        2
Batch1_ATGTGTGCAATAGCGG-1          yes              3                        2
Batch1_CAGTCCTAGCTTTGGT-1          yes              3                        2
Batch1_CGATTGATCTATCCCG-1          yes              3                        2
                          No.Prior.Radiations MRPplus.Overall.Best.Response
Batch1_ACGCAGCAGTCATCCA-1                   2      Progression Disease (PD)
Batch1_ACTGAGTGTTCCGGCA-1                   2      Progression Disease (PD)
Batch1_AGCATACGTAGCTAAA-1                   2      Progression Disease (PD)
Batch1_ATGTGTGCAATAGCGG-1                   2      Progression Disease (PD)
Batch1_CAGTCCTAGCTTTGGT-1                   2      Progression Disease (PD)
Batch1_CGATTGATCTATCCCG-1                   2      Progression Disease (PD)
                          binary_response CD3_score CD3_high Collection.Date
Batch1_ACGCAGCAGTCATCCA-1              PD         1    FALSE            <NA>
Batch1_ACTGAGTGTTCCGGCA-1              PD         1    FALSE            <NA>
Batch1_AGCATACGTAGCTAAA-1              PD         1    FALSE            <NA>
Batch1_ATGTGTGCAATAGCGG-1              PD         1    FALSE            <NA>
Batch1_CAGTCCTAGCTTTGGT-1              PD         1    FALSE            <NA>
Batch1_CGATTGATCTATCCCG-1              PD         1    FALSE            <NA>
                          EGFR.Amplification EGFRvIII..Exon.2.7.deletion
Batch1_ACGCAGCAGTCATCCA-1               <NA>                        <NA>
Batch1_ACTGAGTGTTCCGGCA-1               <NA>                        <NA>
Batch1_AGCATACGTAGCTAAA-1               <NA>                        <NA>
Batch1_ATGTGTGCAATAGCGG-1               <NA>                        <NA>
Batch1_CAGTCCTAGCTTTGGT-1               <NA>                        <NA>
Batch1_CGATTGATCTATCCCG-1               <NA>                        <NA>
                          EGFR.Missense TP53.Frameshift TP53.Missense
Batch1_ACGCAGCAGTCATCCA-1          <NA>            <NA>          <NA>
Batch1_ACTGAGTGTTCCGGCA-1          <NA>            <NA>          <NA>
Batch1_AGCATACGTAGCTAAA-1          <NA>            <NA>          <NA>
Batch1_ATGTGTGCAATAGCGG-1          <NA>            <NA>          <NA>
Batch1_CAGTCCTAGCTTTGGT-1          <NA>            <NA>          <NA>
Batch1_CGATTGATCTATCCCG-1          <NA>            <NA>          <NA>
                          TP53.Nonsense TERT.Promoter.Mutation IDH1.Missense
Batch1_ACGCAGCAGTCATCCA-1          <NA>                   <NA>          <NA>
Batch1_ACTGAGTGTTCCGGCA-1          <NA>                   <NA>          <NA>
Batch1_AGCATACGTAGCTAAA-1          <NA>                   <NA>          <NA>
Batch1_ATGTGTGCAATAGCGG-1          <NA>                   <NA>          <NA>
Batch1_CAGTCCTAGCTTTGGT-1          <NA>                   <NA>          <NA>
Batch1_CGATTGATCTATCCCG-1          <NA>                   <NA>          <NA>
                          IDH2.Missense PTEN.Codon.Deletion PTEN.Frameshift
Batch1_ACGCAGCAGTCATCCA-1          <NA>                <NA>            <NA>
Batch1_ACTGAGTGTTCCGGCA-1          <NA>                <NA>            <NA>
Batch1_AGCATACGTAGCTAAA-1          <NA>                <NA>            <NA>
Batch1_ATGTGTGCAATAGCGG-1          <NA>                <NA>            <NA>
Batch1_CAGTCCTAGCTTTGGT-1          <NA>                <NA>            <NA>
Batch1_CGATTGATCTATCCCG-1          <NA>                <NA>            <NA>
                          PTEN.Interference.of.splice.acceptor.site.in.Intron.6
Batch1_ACGCAGCAGTCATCCA-1                                                  <NA>
Batch1_ACTGAGTGTTCCGGCA-1                                                  <NA>
Batch1_AGCATACGTAGCTAAA-1                                                  <NA>
Batch1_ATGTGTGCAATAGCGG-1                                                  <NA>
Batch1_CAGTCCTAGCTTTGGT-1                                                  <NA>
Batch1_CGATTGATCTATCCCG-1                                                  <NA>
                          CD3_high_low cluster celltype myeloid_cluster
Batch1_ACGCAGCAGTCATCCA-1          Low       0   Lymph1            <NA>
Batch1_ACTGAGTGTTCCGGCA-1          Low       0   Lymph1            <NA>
Batch1_AGCATACGTAGCTAAA-1          Low       0   Lymph1            <NA>
Batch1_ATGTGTGCAATAGCGG-1          Low       0   Lymph1            <NA>
Batch1_CAGTCCTAGCTTTGGT-1          Low       8   Lymph2            <NA>
Batch1_CGATTGATCTATCCCG-1          Low       0   Lymph1            <NA>
                          celltype2 nCount_integrated_sct
Batch1_ACGCAGCAGTCATCCA-1      Teff                     0
Batch1_ACTGAGTGTTCCGGCA-1      Teff                     0
Batch1_AGCATACGTAGCTAAA-1      Teff                     0
Batch1_ATGTGTGCAATAGCGG-1      Teff                     0
Batch1_CAGTCCTAGCTTTGGT-1        NK                     0
Batch1_CGATTGATCTATCCCG-1      Teff                     0
                          nFeature_integrated_sct TGFbELISA TGFbELISA2
Batch1_ACGCAGCAGTCATCCA-1                       0        ND         ND
Batch1_ACTGAGTGTTCCGGCA-1                       0        ND         ND
Batch1_AGCATACGTAGCTAAA-1                       0        ND         ND
Batch1_ATGTGTGCAATAGCGG-1                       0        ND         ND
Batch1_CAGTCCTAGCTTTGGT-1                       0        ND         ND
Batch1_CGATTGATCTATCCCG-1                       0        ND         ND

table(immune_robyn$seurat_clusters)

   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
4642 4509 3785 3260 2549 1569    0 1095  822  608  562  451  368  344  315  291 
  16   17   18 
 142  125    0 

immune_robyn$seurat_clusters <- as.character(immune_robyn$seurat_clusters)

immune_seurat$skull_cluster <- mapvalues(x = colnames(immune_seurat),
                                         from = colnames(immune_robyn),
                                         to = immune_robyn$seurat_clusters)

head(table(immune_seurat$skull_cluster), n = 20)

                      0                       1                      10 
                   4540                    4497                     560 
                     11                      12                      13 
                    443                     360                     340 
                     14                      15                      16 
                    313                     290                     135 
                     17                       2                       3 
                    118                    3720                    3195 
                      4 42-1_AAGGAGCAGCACACAG-1 42-1_AAGGTTCGTAGAGGAA-1 
                   2513                       1                       1 
42-1_ACTGAGTCAATGGACG-1 42-1_AGCGTCGCAAGTTAAG-1 42-1_AGGTCATTCCTAAGTG-1 
                      1                       1                       1 
42-1_ATCCGAAGTCATGCAT-1 42-1_ATGCGATAGTACGCGA-1 
                      1                       1 

table(immune_seurat$skull_cluster)["0"]

   0 
4540 

immune_seurat$skull_cluster <- ifelse(colnames(immune_seurat) %in% colnames(immune_robyn), immune_seurat$skull_cluster, NA)

# Plotting cell type annotations
DimPlot(immune_seurat,
        group.by = "celltype",
        cols = immune_fibro_celltype_col,
        reduction = reduction,
        label = T,
        label.box = T,
        label.size = 3,
        repel = T,
        raster = T,
        raster.dpi = c(1024, 1024),
        pt.size = 3) +
  ggtitle("") +
  theme_classic() +
  NoLegend() +
  NoAxes() +
  coord_fixed(1)

Past versions of unnamed-chunk-2-2.png

Version	Author	Date
9015511	heinin	2024-08-26

DimPlot(immune_seurat,
        group.by = "skull_cluster",
        #cols = immune_fibro_celltype_col,
        reduction = reduction,
        label = T,
        label.box = T,
        label.size = 3,
        repel = T,
        raster = T,
        raster.dpi = c(1024, 1024),
        pt.size = 3) +
  ggtitle("") +
  theme_classic() +
  NoLegend() +
  NoAxes() +
  coord_fixed(1)

DimPlot(tumor_seurat,
        group.by = "celltype",
        cols = tumor_celltype_col,
        reduction = "umap",
        label = T,
        label.box = T,
        label.size = 3,
        repel = T,
        raster = T,
        raster.dpi = c(1024, 1024),
        pt.size = 3) +
  ggtitle("") +
  theme_classic() +
  NoLegend() +
  NoAxes() +
  coord_fixed(1)

rownames(tumor_seurat)[grepl("^TGF", rownames(tumor_seurat))]

 [1] "TGFBR3"   "TGFB2"    "TGFA"     "TGFBRAP1" "TGFBR2"   "TGFBI"   
 [7] "TGFBR1"   "TGFB3"    "TGFB1I1"  "TGFBR3L"  "TGFB1"   

VlnPlot(tumor_seurat, group.by = "TGFbELISA2",
        features = rownames(tumor_seurat)[grepl("^TGF", rownames(tumor_seurat))],
        ncol = 4,
        slot = "counts")

Past versions of unnamed-chunk-3-1.png

Version	Author	Date
9015511	heinin	2024-08-26

# Correlation between measurements and mRNA level
tgfb_elisa <- read.table("/home/hnatri/CART_TGFb/data/ELISA1_TGFb.tsv", sep = "\t", header = T)
tumor_seurat$log10TGFb <- mapvalues(x = tumor_seurat$UPN,
                                    from = tgfb_elisa$UPN,
                                    to = tgfb_elisa$log10TGFb)

layer_data_RNA <- LayerData(tumor_seurat,
                            assay = "RNA",
                            layer = "counts")
layer_data_RNA <- as.data.frame(t(as.matrix(layer_data_RNA)))

layer_data_RNA[1:10, 1:10]

                          OR4F5 OR4F29 OR4F16 SAMD11 NOC2L KLHL17 PLEKHN1 PERM1
Batch1_ACGGGTCTCCAGAGGA-1     0      0      0      0     0      0       0     0
Batch1_ACTTTCAGTCCGACGT-1     0      0      0      0     0      0       0     0
Batch1_AGGGTGACACCTGGTG-1     0      0      0      0     0      0       0     0
Batch1_CACACAAGTACCGTAT-1     0      0      0      0     0      0       0     0
Batch1_CCTTCCCCAAGAAGAG-1     0      0      0      0     1      0       0     0
Batch1_CGATGTACAGGAATCG-1     0      0      0      0     0      0       0     0
Batch1_GAACCTACACCGATAT-1     0      0      0      0     0      0       0     0
Batch1_GAATGAAGTCACCTAA-1     0      0      0      0     2      1       0     0
Batch1_GCAGCCAAGTGTCCCG-1     0      0      0      0     0      0       0     0
Batch1_GTCTTCGAGTCGTTTG-1     0      0      0      0     1      0       0     0
                          HES4 ISG15
Batch1_ACGGGTCTCCAGAGGA-1    0     0
Batch1_ACTTTCAGTCCGACGT-1    0     0
Batch1_AGGGTGACACCTGGTG-1    0     0
Batch1_CACACAAGTACCGTAT-1    0     0
Batch1_CCTTCCCCAAGAAGAG-1    1     0
Batch1_CGATGTACAGGAATCG-1    0     0
Batch1_GAACCTACACCGATAT-1    0     0
Batch1_GAATGAAGTCACCTAA-1    5     0
Batch1_GCAGCCAAGTGTCCCG-1    0     0
Batch1_GTCTTCGAGTCGTTTG-1    1     0

gene_list <- c(rownames(tumor_seurat)[grepl("^TGF", rownames(tumor_seurat))],
               "SMAD2", "SMAD3", "SMAD4")
layer_data_RNA <- layer_data_RNA[, gene_list]

identical(rownames(layer_data_RNA), colnames(tumor_seurat))

[1] TRUE

layer_data_RNA$log10TGFb <- tumor_seurat$log10TGFb
layer_data_RNA$UPN <- tumor_seurat$UPN
layer_data_RNA$TGFbELISA2 <- tumor_seurat$TGFbELISA2
layer_data_RNA_filtered <- layer_data_RNA %>% filter(UPN %in% tgfb_elisa$UPN)

head(layer_data_RNA_filtered)

                          TGFBR3 TGFB2 TGFA TGFBRAP1 TGFBR2 TGFBI TGFBR1 TGFB3
Batch1_AACTTTCTCACCACCT-1      0     0    0        0      0     0      0     0
Batch1_AGATTGCGTAATTGGA-1      0     0    0        0      0     0      0     0
Batch1_AGCTCTCGTATGAATG-1      0     0    0        0      0     0      1     0
Batch1_ATTGGACGTAGTACCT-1      0     0    0        0      0     0      0     0
Batch1_CACCAGGCAATAGCAA-1      0     0    0        0      0     0      0     0
Batch1_CATCAGAGTCCAGTGC-1      0     0    0        0      0     0      0     0
                          TGFB1I1 TGFBR3L TGFB1 SMAD2 SMAD3 SMAD4   log10TGFb
Batch1_AACTTTCTCACCACCT-1       0       0     0     0     0     0 2.072273095
Batch1_AGATTGCGTAATTGGA-1       0       0     0     0     0     0 2.072273095
Batch1_AGCTCTCGTATGAATG-1       0       0     0     0     0     0 2.072273095
Batch1_ATTGGACGTAGTACCT-1       0       0     0     0     0     1 2.072273095
Batch1_CACCAGGCAATAGCAA-1       0       0     0     0     1     0 2.072273095
Batch1_CATCAGAGTCCAGTGC-1       0       0     1     1     0     0 2.072273095
                          UPN TGFbELISA2
Batch1_AACTTTCTCACCACCT-1 185    LowTGFb
Batch1_AGATTGCGTAATTGGA-1 185    LowTGFb
Batch1_AGCTCTCGTATGAATG-1 185    LowTGFb
Batch1_ATTGGACGTAGTACCT-1 185    LowTGFb
Batch1_CACCAGGCAATAGCAA-1 185    LowTGFb
Batch1_CATCAGAGTCCAGTGC-1 185    LowTGFb

layer_data_RNA_filtered %>% ggplot(aes(x = as.numeric(log10TGFb), y = TGFB2, color = TGFbELISA2)) +
  geom_point() +
  theme_bw()

layer_data_RNA_filtered %>% ggplot(aes(x = as.numeric(log10TGFb), y = log2(TGFB2), color = TGFbELISA2)) +
  geom_point() +
  theme_bw()

layer_data_RNA_filtered %>% ggplot(aes(x = as.numeric(log10TGFb), y = TGFBI, color = TGFbELISA2)) +
  geom_point() +
  theme_bw()

layer_data_RNA_filtered %>% ggplot(aes(x = as.numeric(log10TGFb), y = log2(TGFBI), color = TGFbELISA2)) +
  geom_point() +
  theme_bw()

Cell type markers for the immune subset

DefaultAssay(immune_seurat) <- "RNA"

# Dropping MT and RP genes before calling markers
RBMTgenes <- grep(pattern = "^RP[SL]|^MRP[SL]|^MT-",
                  x = rownames(immune_seurat@assays$RNA@data),
                  value = TRUE, invert = TRUE)
immune_seurat <- subset(immune_seurat, features = RBMTgenes)

# Top markers for each cluster
markers <- presto::wilcoxauc(immune_seurat,
                             group_by = "celltype",
                             assay = "data",
                             seurat_assay = "RNA")

top_markers <- markers %>% group_by(group) %>% slice_max(order_by = auc, n = 2)

FeaturePlot(immune_seurat,
            features = top_markers$feature,
            ncol = 6,
            reduction = "umap",
            raster = T,
            cols = c("gray89", "tomato3")) &
  coord_fixed(ratio = 1) &
  theme_bw() &
  NoLegend() &
  manuscript_theme

Past versions of unnamed-chunk-4-1.png

Version	Author	Date
9015511	heinin	2024-08-26

top_markers <- markers %>%  group_by(group) %>% slice_max(order_by = auc, n = 5)

# seurat_object, plot_features, group_var, group_colors, column_title, km=5, row.order = NULL
create_dotplot_heatmap_horizontal(seurat_object = immune_seurat,
                                  plot_features = unique(top_markers$feature),
                                  group_var = "celltype",
                                  group_colors = immune_fibro_celltype_col,
                                  column_title = "",
                                  km = 5,
                                  col.order = NULL)

Past versions of unnamed-chunk-5-1.png

Version	Author	Date
9015511	heinin	2024-08-26

Past versions of unnamed-chunk-5-2.png

Version	Author	Date
9015511	heinin	2024-08-26

Differential gene expression

unique(immune_seurat$TGFbELISA2)

[1] "ND"       "LowTGFb"  "HighTGFb"

# For each cluster, top DEGs between TGFb high and low
DEG_TGFb <- lapply(unique(immune_seurat$celltype), function(xx){
  data_subset <- subset(immune_seurat, subset = celltype == xx)
  Idents(data_subset) <- data_subset$TGFbELISA2
  if (all((c("HighTGFb", "LowTGFb") %in% data_subset$TGFbELISA2) == c(T, T))){
    markers <- FindMarkers(data_subset,
                           ident.1 = "HighTGFb",
                           ident.2 = "LowTGFb",
                           assay = "RNA",
                           verbose = F)
    markers$feature <- rownames(markers)
    markers$celltype <- xx
    
    return(markers)
  } else {
    return(NULL)
  }
})
names(DEG_TGFb) <- unique(immune_seurat$celltype)
DEG_TGFb[sapply(DEG_TGFb, is.null)] <- NULL

DEG_TGFb_df <- as.data.frame(do.call(rbind, DEG_TGFb))

# Distribution of p-values and log2FC
hist(DEG_TGFb_df$p_val)

Past versions of unnamed-chunk-6-1.png

Version	Author	Date
9015511	heinin	2024-08-26

hist(DEG_TGFb_df$avg_log2FC)

Past versions of unnamed-chunk-6-2.png

Version	Author	Date
9015511	heinin	2024-08-26

DEG_TGFb_df_sig <- DEG_TGFb_df %>%
  filter(p_val < 0.01,
         abs(avg_log2FC) > 2,
         (pct.1 > 0.50 | pct.2 > 0.50))

dim(DEG_TGFb_df_sig)

[1] 235   7

# Saving to a file
write.table(DEG_TGFb_df_sig,
            "/scratch/hnatri/CART/TGFb_high_vs_low_DEGs_immune_sig.tsv",
            sep = "\t", quote = F, row.names = F)

# Plotting
table(DEG_TGFb_df_sig$celltype) %>% as.data.frame() %>%
  ggplot(aes(x = reorder(Var1, -Freq), y = Freq, fill = Var1)) +
    geom_bar(stat = "identity") +
    scale_fill_manual(values = immune_fibro_celltype_col) +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
    NoLegend() +
    xlab("Cell type") +
    ylab("# DEGs, TGFb high vs. low")

Pathways enrichment

# Adding HALLMARK pathways from another object
immune_fibro_hallmark <- readRDS("/tgen_labs/banovich/BCTCSF/Heini/13384_Tumor/13384_immune_fibro_GSEA_C2_H.rds")
colnames(immune_fibro_hallmark@meta.data)[grep("HALLMARK", colnames(immune_fibro_hallmark@meta.data))]

 [1] "HALLMARK_ADIPOGENESIS"                     
 [2] "HALLMARK_ALLOGRAFT_REJECTION"              
 [3] "HALLMARK_ANDROGEN_RESPONSE"                
 [4] "HALLMARK_ANGIOGENESIS"                     
 [5] "HALLMARK_APICAL_JUNCTION"                  
 [6] "HALLMARK_APICAL_SURFACE"                   
 [7] "HALLMARK_APOPTOSIS"                        
 [8] "HALLMARK_BILE_ACID_METABOLISM"             
 [9] "HALLMARK_CHOLESTEROL_HOMEOSTASIS"          
[10] "HALLMARK_COAGULATION"                      
[11] "HALLMARK_COMPLEMENT"                       
[12] "HALLMARK_DNA_REPAIR"                       
[13] "HALLMARK_E2F_TARGETS"                      
[14] "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION"
[15] "HALLMARK_ESTROGEN_RESPONSE_EARLY"          
[16] "HALLMARK_ESTROGEN_RESPONSE_LATE"           
[17] "HALLMARK_FATTY_ACID_METABOLISM"            
[18] "HALLMARK_G2M_CHECKPOINT"                   
[19] "HALLMARK_GLYCOLYSIS"                       
[20] "HALLMARK_HEDGEHOG_SIGNALING"               
[21] "HALLMARK_HEME_METABOLISM"                  
[22] "HALLMARK_HYPOXIA"                          
[23] "HALLMARK_IL2_STAT5_SIGNALING"              
[24] "HALLMARK_IL6_JAK_STAT3_SIGNALING"          
[25] "HALLMARK_INFLAMMATORY_RESPONSE"            
[26] "HALLMARK_INTERFERON_ALPHA_RESPONSE"        
[27] "HALLMARK_INTERFERON_GAMMA_RESPONSE"        
[28] "HALLMARK_KRAS_SIGNALING_DN"                
[29] "HALLMARK_KRAS_SIGNALING_UP"                
[30] "HALLMARK_MITOTIC_SPINDLE"                  
[31] "HALLMARK_MTORC1_SIGNALING"                 
[32] "HALLMARK_MYC_TARGETS_V1"                   
[33] "HALLMARK_MYC_TARGETS_V2"                   
[34] "HALLMARK_MYOGENESIS"                       
[35] "HALLMARK_NOTCH_SIGNALING"                  
[36] "HALLMARK_OXIDATIVE_PHOSPHORYLATION"        
[37] "HALLMARK_P53_PATHWAY"                      
[38] "HALLMARK_PANCREAS_BETA_CELLS"              
[39] "HALLMARK_PEROXISOME"                       
[40] "HALLMARK_PI3K_AKT_MTOR_SIGNALING"          
[41] "HALLMARK_PROTEIN_SECRETION"                
[42] "HALLMARK_REACTIVE_OXYGEN_SPECIES_PATHWAY"  
[43] "HALLMARK_SPERMATOGENESIS"                  
[44] "HALLMARK_TGF_BETA_SIGNALING"               
[45] "HALLMARK_TNFA_SIGNALING_VIA_NFKB"          
[46] "HALLMARK_UNFOLDED_PROTEIN_RESPONSE"        
[47] "HALLMARK_UV_RESPONSE_DN"                   
[48] "HALLMARK_UV_RESPONSE_UP"                   
[49] "HALLMARK_WNT_BETA_CATENIN_SIGNALING"       
[50] "HALLMARK_XENOBIOTIC_METABOLISM"            

for(i in colnames(immune_fibro_hallmark@meta.data)[grep("HALLMARK", colnames(immune_fibro_hallmark@meta.data))]){
  immune_seurat@meta.data[,i] <- mapvalues(x = rownames(immune_seurat@meta.data),
                                           from = rownames(immune_fibro_hallmark@meta.data),
                                           to = immune_fibro_hallmark@meta.data[,i])
}

rm(immune_fibro_hallmark)

tumors_hallmark <- readRDS("/tgen_labs/banovich/BCTCSF/Heini/13384_Tumor/13384_tumors_GSEA_C2_H.rds")
colnames(tumors_hallmark@meta.data)[grep("HALLMARK", colnames(tumors_hallmark@meta.data))]

 [1] "HALLMARK_ADIPOGENESIS"                     
 [2] "HALLMARK_ALLOGRAFT_REJECTION"              
 [3] "HALLMARK_ANDROGEN_RESPONSE"                
 [4] "HALLMARK_ANGIOGENESIS"                     
 [5] "HALLMARK_APICAL_JUNCTION"                  
 [6] "HALLMARK_APICAL_SURFACE"                   
 [7] "HALLMARK_APOPTOSIS"                        
 [8] "HALLMARK_BILE_ACID_METABOLISM"             
 [9] "HALLMARK_CHOLESTEROL_HOMEOSTASIS"          
[10] "HALLMARK_COAGULATION"                      
[11] "HALLMARK_COMPLEMENT"                       
[12] "HALLMARK_DNA_REPAIR"                       
[13] "HALLMARK_E2F_TARGETS"                      
[14] "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION"
[15] "HALLMARK_ESTROGEN_RESPONSE_EARLY"          
[16] "HALLMARK_ESTROGEN_RESPONSE_LATE"           
[17] "HALLMARK_FATTY_ACID_METABOLISM"            
[18] "HALLMARK_G2M_CHECKPOINT"                   
[19] "HALLMARK_GLYCOLYSIS"                       
[20] "HALLMARK_HEDGEHOG_SIGNALING"               
[21] "HALLMARK_HEME_METABOLISM"                  
[22] "HALLMARK_HYPOXIA"                          
[23] "HALLMARK_IL2_STAT5_SIGNALING"              
[24] "HALLMARK_IL6_JAK_STAT3_SIGNALING"          
[25] "HALLMARK_INFLAMMATORY_RESPONSE"            
[26] "HALLMARK_INTERFERON_ALPHA_RESPONSE"        
[27] "HALLMARK_INTERFERON_GAMMA_RESPONSE"        
[28] "HALLMARK_KRAS_SIGNALING_DN"                
[29] "HALLMARK_KRAS_SIGNALING_UP"                
[30] "HALLMARK_MITOTIC_SPINDLE"                  
[31] "HALLMARK_MTORC1_SIGNALING"                 
[32] "HALLMARK_MYC_TARGETS_V1"                   
[33] "HALLMARK_MYC_TARGETS_V2"                   
[34] "HALLMARK_MYOGENESIS"                       
[35] "HALLMARK_NOTCH_SIGNALING"                  
[36] "HALLMARK_OXIDATIVE_PHOSPHORYLATION"        
[37] "HALLMARK_P53_PATHWAY"                      
[38] "HALLMARK_PANCREAS_BETA_CELLS"              
[39] "HALLMARK_PEROXISOME"                       
[40] "HALLMARK_PI3K_AKT_MTOR_SIGNALING"          
[41] "HALLMARK_PROTEIN_SECRETION"                
[42] "HALLMARK_REACTIVE_OXYGEN_SPECIES_PATHWAY"  
[43] "HALLMARK_SPERMATOGENESIS"                  
[44] "HALLMARK_TGF_BETA_SIGNALING"               
[45] "HALLMARK_TNFA_SIGNALING_VIA_NFKB"          
[46] "HALLMARK_UNFOLDED_PROTEIN_RESPONSE"        
[47] "HALLMARK_UV_RESPONSE_DN"                   
[48] "HALLMARK_UV_RESPONSE_UP"                   
[49] "HALLMARK_WNT_BETA_CATENIN_SIGNALING"       
[50] "HALLMARK_XENOBIOTIC_METABOLISM"            

for(i in colnames(tumors_hallmark@meta.data)[grep("HALLMARK", colnames(tumors_hallmark@meta.data))]){
  tumor_seurat@meta.data[,i] <- mapvalues(x = rownames(tumor_seurat@meta.data),
                                          from = rownames(tumors_hallmark@meta.data),
                                          to = tumors_hallmark@meta.data[,i])
}

rm(tumors_hallmark)

#saveRDS(immune_seurat, "/tgen_labs/banovich/BCTCSF/Heini/13384_Tumor/13384_immune_GSEA_C2_H.rds")
#saveRDS(tumor_seurat, "/tgen_labs/banovich/BCTCSF/Heini/13384_Tumor/13384_tumor_GSEA_C2_H.rds")
immune_seurat <- readRDS("/tgen_labs/banovich/BCTCSF/Heini/13384_Tumor/13384_immune_GSEA_C2_H.rds")
tumor_seurat <- readRDS("/tgen_labs/banovich/BCTCSF/Heini/13384_Tumor/13384_tumor_GSEA_C2_H.rds")

colnames(immune_seurat@meta.data)[grep("HALLMARK", colnames(immune_seurat@meta.data))]

 [1] "HALLMARK_ADIPOGENESIS"                     
 [2] "HALLMARK_ALLOGRAFT_REJECTION"              
 [3] "HALLMARK_ANDROGEN_RESPONSE"                
 [4] "HALLMARK_ANGIOGENESIS"                     
 [5] "HALLMARK_APICAL_JUNCTION"                  
 [6] "HALLMARK_APICAL_SURFACE"                   
 [7] "HALLMARK_APOPTOSIS"                        
 [8] "HALLMARK_BILE_ACID_METABOLISM"             
 [9] "HALLMARK_CHOLESTEROL_HOMEOSTASIS"          
[10] "HALLMARK_COAGULATION"                      
[11] "HALLMARK_COMPLEMENT"                       
[12] "HALLMARK_DNA_REPAIR"                       
[13] "HALLMARK_E2F_TARGETS"                      
[14] "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION"
[15] "HALLMARK_ESTROGEN_RESPONSE_EARLY"          
[16] "HALLMARK_ESTROGEN_RESPONSE_LATE"           
[17] "HALLMARK_FATTY_ACID_METABOLISM"            
[18] "HALLMARK_G2M_CHECKPOINT"                   
[19] "HALLMARK_GLYCOLYSIS"                       
[20] "HALLMARK_HEDGEHOG_SIGNALING"               
[21] "HALLMARK_HEME_METABOLISM"                  
[22] "HALLMARK_HYPOXIA"                          
[23] "HALLMARK_IL2_STAT5_SIGNALING"              
[24] "HALLMARK_IL6_JAK_STAT3_SIGNALING"          
[25] "HALLMARK_INFLAMMATORY_RESPONSE"            
[26] "HALLMARK_INTERFERON_ALPHA_RESPONSE"        
[27] "HALLMARK_INTERFERON_GAMMA_RESPONSE"        
[28] "HALLMARK_KRAS_SIGNALING_DN"                
[29] "HALLMARK_KRAS_SIGNALING_UP"                
[30] "HALLMARK_MITOTIC_SPINDLE"                  
[31] "HALLMARK_MTORC1_SIGNALING"                 
[32] "HALLMARK_MYC_TARGETS_V1"                   
[33] "HALLMARK_MYC_TARGETS_V2"                   
[34] "HALLMARK_MYOGENESIS"                       
[35] "HALLMARK_NOTCH_SIGNALING"                  
[36] "HALLMARK_OXIDATIVE_PHOSPHORYLATION"        
[37] "HALLMARK_P53_PATHWAY"                      
[38] "HALLMARK_PANCREAS_BETA_CELLS"              
[39] "HALLMARK_PEROXISOME"                       
[40] "HALLMARK_PI3K_AKT_MTOR_SIGNALING"          
[41] "HALLMARK_PROTEIN_SECRETION"                
[42] "HALLMARK_REACTIVE_OXYGEN_SPECIES_PATHWAY"  
[43] "HALLMARK_SPERMATOGENESIS"                  
[44] "HALLMARK_TGF_BETA_SIGNALING"               
[45] "HALLMARK_TNFA_SIGNALING_VIA_NFKB"          
[46] "HALLMARK_UNFOLDED_PROTEIN_RESPONSE"        
[47] "HALLMARK_UV_RESPONSE_DN"                   
[48] "HALLMARK_UV_RESPONSE_UP"                   
[49] "HALLMARK_WNT_BETA_CATENIN_SIGNALING"       
[50] "HALLMARK_XENOBIOTIC_METABOLISM"            

# Pathway names
pathway_names <- c(colnames(immune_seurat@meta.data)[grep("REACTOME", colnames(immune_seurat@meta.data))],
                   colnames(immune_seurat@meta.data)[grep("KEGG", colnames(immune_seurat@meta.data))],
                   colnames(immune_seurat@meta.data)[grep("BIOCARTA", colnames(immune_seurat@meta.data))],
                   colnames(immune_seurat@meta.data)[grep("HALLMARK", colnames(immune_seurat@meta.data))])

table(sapply(strsplit(pathway_names, split='_', fixed=TRUE), `[`, 1))

BIOCARTA HALLMARK     KEGG REACTOME 
     291       50      186     1602 

# Plotting results for all immune cells
gsea_res <- immune_seurat@meta.data %>%
  dplyr::select(all_of(c("TGFbELISA2", pathway_names)))

gsea_res <- gsea_res %>% filter(TGFbELISA2 %in% c("HighTGFb", "LowTGFb"))
output_immune <- data.frame(getSignificance(gsea_res, group = "TGFbELISA2", fit = "Wilcoxon"))
output_immune$pathways <- rownames(output_immune)
output_immune <- output_immune %>% filter(FDR < 0.01) %>% arrange(FDR)

VlnPlot(immune_seurat,
        group.by = "TGFbELISA2",
        features = c("REACTOME_SIGNALING_BY_TGFB_FAMILY_MEMBERS"),
        pt.size = 0)

Plotting TGFb pathways in all immune cells

# Plotting TGFb pathways
output_immune_tgf <- output_immune %>% filter(pathways %in% output_immune$pathways[grep("TGF", output_immune$pathways)])

output_immune_tgf %>%
  mutate(delta = median.HighTGFb - median.LowTGFb,
         sign = sign(delta),
         signstr = if_else(sign == 1, "HighTGFb", "LowTGFb")) %>%
  ggplot(aes(x = delta, y = reorder(pathways, delta), fill = signstr)) +
    geom_bar(stat = "identity") +
    #geom_col(width = 0.85) +
    scale_fill_manual(values = c("orangered1", "royalblue3")) +
    theme_classic() +
    manuscript_theme +
    ylab("") +
    xlab(expression(Delta ~ "median enrichment score"))

    #coord_flip()

# Selecting a subset of pathways for plotting
output_immune_plot <- output_immune %>% head(n=20)

output_immune_plot %>%
  mutate(delta = median.HighTGFb - median.LowTGFb,
         sign = sign(delta),
         signstr = if_else(sign == 1, "HighTGFb", "LowTGFb")) %>%
  ggplot(aes(x = delta, y = reorder(pathways, delta), fill = signstr)) +
    geom_bar(stat = "identity") +
    #geom_col(width = 0.85) +
    scale_fill_manual(values = c("orangered1", "royalblue3")) +
    theme_classic() +
    manuscript_theme +
    ylab("") +
    xlab(expression(Delta ~ "median enrichment score"))

    #coord_flip()

write.table(output_immune, "/scratch/hnatri/CART/TGFb_high_vs_low_GSEA_sig.tsv",
            quote = F, row.names = F, sep = "\t")

Myeloid only GSEA

# Myeloid celltypes only
myeloid <- subset(immune_seurat, subset = celltype %in% c(paste0("M", seq(1, 9)), "B1", "N1"))
gsea_res_myeloid <- myeloid@meta.data %>%
  dplyr::select(c("orig.ident", "TGFbELISA2", pathway_names))

gsea_res_myeloid <- gsea_res_myeloid %>% filter(TGFbELISA2 %in% c("HighTGFb", "LowTGFb"))
output_myeloid <- data.frame(getSignificance(gsea_res_myeloid, group = "TGFbELISA2", fit = "Wilcoxon"))
output_myeloid$pathways <- rownames(output_myeloid)
output_myeloid <- output_myeloid %>% filter(FDR < 0.01) %>% arrange(FDR)

write.table(output_myeloid, "/scratch/hnatri/CART/TGFb_high_vs_low_myeloid_GSEA_sig.tsv",
            quote = F, row.names = F, sep = "\t")

output_myeloid$pathways[grep("INTERFERON", output_myeloid$pathways)]

[1] "REACTOME_INTERFERON_GAMMA_SIGNALING"                             
[2] "REACTOME_INTERFERON_SIGNALING"                                   
[3] "REACTOME_INTERFERON_ALPHA_BETA_SIGNALING"                        
[4] "REACTOME_DDX58_IFIH1_MEDIATED_INDUCTION_OF_INTERFERON_ALPHA_BETA"

TGFbhigh_enrichment <- myeloid@meta.data %>% filter(TGFbELISA2 == "HighTGFb") %>%
  dplyr::select(REACTOME_INTERFERON_GAMMA_SIGNALING) %>% unlist() %>% as.numeric()
TGFblow_enrichment <- myeloid@meta.data %>% filter(TGFbELISA2 == "LowTGFb") %>%
  dplyr::select(REACTOME_INTERFERON_GAMMA_SIGNALING) %>% unlist() %>% as.numeric()

summary(TGFbhigh_enrichment)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   2365    4199    4597    4602    4987    6855 

summary(TGFblow_enrichment)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3069    4620    5069    5062    5516    7241 

t.test(TGFbhigh_enrichment, TGFblow_enrichment)

    Welch Two Sample t-test

data:  TGFbhigh_enrichment and TGFblow_enrichment
t = -28.774, df = 6107, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -490.8579 -428.2405
sample estimates:
mean of x mean of y 
 4602.029  5061.579 

output_myeloid[which(output_myeloid$pathways == "REACTOME_INTERLEUKIN_27_SIGNALING"),]

   W.statistic      p.value          FDR median.LowTGFb median.HighTGFb
31     3189886 6.53562e-101 1.339149e-97       1370.491        2315.554
                            pathways
31 REACTOME_INTERLEUKIN_27_SIGNALING

output_myeloid[which(output_myeloid$pathways == "REACTOME_INTERFERON_GAMMA_SIGNALING"),]

  W.statistic       p.value           FDR median.LowTGFb median.HighTGFb
2     2794825 2.795529e-161 5.809109e-158       4596.855        5069.075
                             pathways
2 REACTOME_INTERFERON_GAMMA_SIGNALING

VlnPlot(myeloid,
        features = "REACTOME_INTERFERON_GAMMA_SIGNALING",
        group.by = "TGFbELISA2",
        pt.size = 0)

# Selecting a subset of pathways for plotting
output_plot_myeloid <- output_myeloid %>% head(n=60)

output_plot_myeloid %>%
  mutate(delta = median.HighTGFb - median.LowTGFb,
         sign = sign(delta),
         signstr = if_else(sign == 1, "HighTGFb", "LowTGFb")) %>%
  ggplot(aes(x = delta, y = reorder(pathways, delta), fill = signstr)) +
    geom_bar(stat = "identity") +
    #geom_col(width = 0.85) +
    scale_fill_manual(values = c("orangered1", "royalblue3")) +
    theme_classic() +
    manuscript_theme +
    ylab("") +
    xlab(expression(Delta ~ "median enrichment score"))

    #coord_flip()

Lymphoid only GSEA

Session information

sessionInfo()

R version 4.3.0 (2023-04-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] circlize_0.4.15             plyr_1.8.8                 
 [3] ComplexHeatmap_2.18.0       viridis_0.6.3              
 [5] viridisLite_0.4.2           RColorBrewer_1.1-3         
 [7] escape_1.12.0               scProportionTest_0.0.0.9000
 [9] patchwork_1.1.2             ggrepel_0.9.3              
[11] lubridate_1.9.2             forcats_1.0.0              
[13] stringr_1.5.0               purrr_1.0.1                
[15] readr_2.1.4                 tidyr_1.3.0                
[17] tibble_3.2.1                ggplot2_3.4.2              
[19] tidyverse_2.0.0             dplyr_1.1.2                
[21] googlesheets4_1.1.0         Seurat_5.0.1               
[23] SeuratObject_5.0.1          sp_1.6-1                   
[25] workflowr_1.7.1            

loaded via a namespace (and not attached):
  [1] fs_1.6.2                    matrixStats_1.0.0          
  [3] GSVA_1.50.5                 spatstat.sparse_3.0-1      
  [5] bitops_1.0-7                doParallel_1.0.17          
  [7] httr_1.4.6                  tools_4.3.0                
  [9] sctransform_0.4.1           backports_1.4.1            
 [11] utf8_1.2.3                  R6_2.5.1                   
 [13] HDF5Array_1.30.1            lazyeval_0.2.2             
 [15] uwot_0.1.14                 GetoptLong_1.0.5           
 [17] rhdf5filters_1.14.1         withr_2.5.0                
 [19] gridExtra_2.3               progressr_0.13.0           
 [21] cli_3.6.1                   Biobase_2.62.0             
 [23] Cairo_1.6-0                 spatstat.explore_3.2-1     
 [25] fastDummies_1.7.3           labeling_0.4.2             
 [27] sass_0.4.6                  spatstat.data_3.0-1        
 [29] ggridges_0.5.4              pbapply_1.7-0              
 [31] parallelly_1.36.0           limma_3.58.1               
 [33] rstudioapi_0.14             RSQLite_2.3.1              
 [35] shape_1.4.6                 generics_0.1.3             
 [37] ica_1.0-3                   spatstat.random_3.1-5      
 [39] Matrix_1.6-5                ggbeeswarm_0.7.2           
 [41] fansi_1.0.4                 S4Vectors_0.40.2           
 [43] abind_1.4-5                 lifecycle_1.0.3            
 [45] whisker_0.4.1               yaml_2.3.7                 
 [47] SummarizedExperiment_1.32.0 rhdf5_2.46.1               
 [49] SparseArray_1.2.3           Rtsne_0.16                 
 [51] blob_1.2.4                  promises_1.2.0.1           
 [53] crayon_1.5.2                miniUI_0.1.1.1             
 [55] lattice_0.21-8              beachmat_2.18.1            
 [57] msigdbr_7.5.1               cowplot_1.1.1              
 [59] annotate_1.80.0             KEGGREST_1.42.0            
 [61] magick_2.7.4                pillar_1.9.0               
 [63] knitr_1.43                  GenomicRanges_1.54.1       
 [65] rjson_0.2.21                future.apply_1.11.0        
 [67] codetools_0.2-19            leiden_0.4.3               
 [69] glue_1.6.2                  getPass_0.2-4              
 [71] data.table_1.14.8           vctrs_0.6.2                
 [73] png_0.1-8                   spam_2.9-1                 
 [75] cellranger_1.1.0            gtable_0.3.3               
 [77] cachem_1.0.8                xfun_0.39                  
 [79] S4Arrays_1.2.0              mime_0.12                  
 [81] survival_3.5-5              gargle_1.4.0               
 [83] SingleCellExperiment_1.24.0 iterators_1.0.14           
 [85] statmod_1.5.0               ellipsis_0.3.2             
 [87] fitdistrplus_1.1-11         ROCR_1.0-11                
 [89] nlme_3.1-162                bit64_4.0.5                
 [91] RcppAnnoy_0.0.20            GenomeInfoDb_1.38.5        
 [93] rprojroot_2.0.3             bslib_0.4.2                
 [95] irlba_2.3.5.1               vipor_0.4.5                
 [97] KernSmooth_2.23-21          colorspace_2.1-0           
 [99] BiocGenerics_0.48.1         DBI_1.1.3                  
[101] ggrastr_1.0.2               UCell_2.6.2                
[103] tidyselect_1.2.0            processx_3.8.1             
[105] curl_5.0.0                  bit_4.0.5                  
[107] compiler_4.3.0              git2r_0.32.0               
[109] graph_1.80.0                BiocNeighbors_1.20.2       
[111] DelayedArray_0.28.0         plotly_4.10.2              
[113] scales_1.2.1                lmtest_0.9-40              
[115] callr_3.7.3                 digest_0.6.31              
[117] goftest_1.2-3               presto_1.0.0               
[119] spatstat.utils_3.0-3        rmarkdown_2.22             
[121] XVector_0.42.0              htmltools_0.5.5            
[123] pkgconfig_2.0.3             sparseMatrixStats_1.14.0   
[125] MatrixGenerics_1.14.0       highr_0.10                 
[127] fastmap_1.1.1               GlobalOptions_0.1.2        
[129] rlang_1.1.1                 htmlwidgets_1.6.2          
[131] shiny_1.7.4                 DelayedMatrixStats_1.24.0  
[133] farver_2.1.1                jquerylib_0.1.4            
[135] zoo_1.8-12                  jsonlite_1.8.5             
[137] BiocParallel_1.36.0         BiocSingular_1.18.0        
[139] RCurl_1.98-1.12             magrittr_2.0.3             
[141] GenomeInfoDbData_1.2.11     dotCall64_1.0-2            
[143] Rhdf5lib_1.24.1             munsell_0.5.0              
[145] Rcpp_1.0.10                 babelgene_22.9             
[147] reticulate_1.29             stringi_1.7.12             
[149] zlibbioc_1.48.0             MASS_7.3-60                
[151] parallel_4.3.0              listenv_0.9.0              
[153] deldir_1.0-9                Biostrings_2.70.1          
[155] splines_4.3.0               tensor_1.5                 
[157] hms_1.1.3                   ps_1.7.5                   
[159] igraph_1.4.3                spatstat.geom_3.2-1        
[161] RcppHNSW_0.5.0              reshape2_1.4.4             
[163] stats4_4.3.0                ScaledMatrix_1.10.0        
[165] XML_3.99-0.14               evaluate_0.21              
[167] foreach_1.5.2               tzdb_0.4.0                 
[169] httpuv_1.6.11               RANN_2.6.1                 
[171] polyclip_1.10-4             clue_0.3-64                
[173] future_1.32.0               scattermore_1.2            
[175] rsvd_1.0.5                  broom_1.0.4                
[177] xtable_1.8-4                RSpectra_0.16-1            
[179] later_1.3.1                 googledrive_2.1.0          
[181] beeswarm_0.4.0              memoise_2.0.1              
[183] AnnotationDbi_1.64.1        IRanges_2.36.0             
[185] cluster_2.1.4               timechange_0.2.0           
[187] globals_0.16.2              GSEABase_1.64.0